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ABSTRACT 


Multi-Grid  Algorithms  with  Applications 
to  Elliptic  Boundary -Value  Problems 

Craig  C.  Douglas 
Yale  University,  1983 

This  work  is  primarily  concerned  with  solving  the  large  sparse  linear  systems  which 
arise  in  connection  with  finite-element  or  finite-difference  procedures  for  solving  self- 
adjoint  elliptic  boundary-value  problems.  These  problems  can  be  expressed  in  terms  of 
abstract  variational  problems  on  Hilbert  spaces.  Our  (multi-grid)  schemes  involve  a 
sequence  of  auxiliary  finite-dimensional  spaces  which  do  not  have  to  be  nested.  We 
approximate  the  solution  using  the  largest  (finite-dimensional)  space.  These  schemes  are 
recursive  in  nature:  they  combine  smoothing  iterations  in  a  space  with  solving  one  or 
more  correction  problems  using  smaller  spaces.  Under  certain  circumstances,  the  solution 
to  a  problem  can  be  approximated  well  using  smaller  spaces.  Since  the  smaller  spaces  are 
required  to  have  geometrically  fewer  unknowns  than  the  largest  space,  the  savings  in 
computation  can  be  substantial.  '"Ip  fact,  we  prove  that  these  procedures  are  optimal 
order  under  appropriate  conditions.^JOur  general  theory  is  discretization  independent  and 
can  be  applied  to  problems  which  do  not  arise  from  partial  differential  equations. 


As  examples,  we  consider  three  particular  discretizations  of  variable  coefficient  self- 
adjoint  second  order  elliptic  boundary-value  problems.  The  first  is  a  finite-element 
discretization  on  a  convex  domain  in  two  dimensions.  The  second  is  a  finite-difference 
discretization  in  one  dimension.  The  last  is  a  finite-difference  discretization  on  the  unit 
square. 
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1:  General  Theory 


In  this  section,  we  discuss  the  approximate  solution  of  an  abstract  elliptic  variational 
problem.  Our  scheme  involves  a  sequence  of  finite-dimensional  spaces  .Mj,  j  **  1,  2,...  ,  k. 
We  approximate  the  solution  using  the  largest  space.  Under  certain  circumstances,  the 
solution  to  a  problem  can  be  approximated  well  using  smaller  spaces.  Since  we  require 
the  smaller  spaces  to  have  geometrically  fewer  unknowns  than  the  largest  one,  the  savings 
in  computation  can  be  substantial.  In  fact,  we  prove  that  these  procedures  are  optimal 
order  under  appropriate  conditions.  While  this  theory  is  applied  to  solving  the  large 
sparse  linear  systems  which  arise  in  connection  with  finite-element  or  finite-difference 
procedures  for  solving  self-adjoint  elliptic  boundary-value  problems  in  the  next  section,  it 
can  also  be  applied  to  problems  which  do  not  arise  from  partial  differential  equations. 

Our  k-level  scheme  is  related  to  the  multi-grid  techniques  used  by  Bank  and 
Dupont  [6],  which  is  related  to  the  techniques  of  Brandt  [10],  Bakhvalov  [5], 
Federenko  [14,  15],  Nicolaides  [25,  26],  and  Hackbusch  [18,  19,  20].  The  earlier  proofs  are 
for  particular  discretizations  of  model  elliptic  boundary- value  problems.  Their  domains 
are  covered  by  meshes  or  triangulations  which  are  refined  uniformly.  Only  Van 
Rosendale's  proof  [29]  allows  nonuniformly  refined  domains.  The  proofs  here  use  abstract 
function  space  arguments  which  make  no  reference  to  the  particular  discretization, 
domain,  or  method  of  refinement.  Further,  we  do  not  require  the  solution  spaces  to  be 
nested  as  in  the  proofs  of  the  cited  references. 

.Assume  we  are  given  a  triple, 

{ H,  a(u,v),  f(v)},  (1.1) 

where  H  is  a  Hilbert  space  with  norm  j|-j|,  a(u,v)  is  a  continuous  symmetric  real  valued 
bilinear  form  on  H  ®  H,  and  f(v):H  — ►  R  is  a  continuous  real  valued  linear  functional. 
F urthermore,  we  assume  that  there  exists  a  constant  aQ  >  0  such  that 

a(v,v)  >  a0j|v||2  for  all  v  6  H. 

The  bilinear  form  a( induces  the  energy  norm 
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IIMII2  35  a(«.4 

We  seek  an  approximation  to  the  solution  of 

Problem  1.1:  Given  {H,  a(u,v),  f(v)},  find  u  6  H  such  that 
a(u,v)  =»  f(v)  for  all  v  6  H. 


Problem  1.1  has  a  unique  solution  (see  Ciarlet  [1 1 J). 

We  now  consider  the  finite-dimensional  approximation  of  Problem  1.1.  Let  M., 
j  >  1,  be  a  sequence  of  N^-dimensional  spaces.  Associated  with  each  space  .Mj  is  a 
continuous,  symmetric,  positive-definite  bilinear  form  aj(u,v)  and  a  continuous,  bounded 
linear  form  f-(v).  We  require  that 

Nj  ~  for  some  a  >  2.  (1.2) 

We  will  see  that  <7  is  important:  when  tr  >  2,  we  can  always  construct  optimal  order 
algorithm!  to  approximate  the  solution  to  Problem  1.1. 

We  assume  that  linear  operators  exist  which  project  H  onto  M.  and  inject  into  H 
for  any  j  >  1: 


P;:  H  X;  and 

onto  j 


V  Mj  r-i>  H- 


(1.3) 


For  j  >  1,  the  linear  operators  defined  by 


Ej:  and 


(1.4) 


Rj:  M j 


interpolate  between  adjacent  solution  spaces.  One  natural  definition  of  Ej  and  Rj  is 


E; 


Pjij-l  “d  Ri  “  Pj-lV 


(15) 


It  is  natural  to  define  R.  as  the  adjoint  of  Ej.- 


Figure  1-1:  Space  Operators 


c 


In  this  case, 

H  -  £-iPj  -  Pj-i'j  -  Ri- 

If  we  assume  that  p-  =*  i- ,  then  all  of  the  operators  can  be  defined  in  terms  of  the 
injection  operators  i-.  However,  we  will  have  occasion  to  use  more  general  operators  than 
those  in  (1.5).  Figure  1-1  shows  the  relationship  between  these  operators  and  the  various 
spaces. 

For  j  >  1,  the  finite-dimensional  approximations  of  Problem  1.1  are 

Problem  1.2:  Given  {.Mj,  a^u.v),  f(v)},  find  Uj  6  .Mj  such  that 
3j(u.,v)  =  f(v)  for  all  v  6  M-}. 

Associated  with  each  space  Mj  are  eigenvalues  and  eigenfunctions  (eigenvectors) 
1  <  i  <  Nj,  satisfying 

aj(T’^^)  ~  M^j  ^or  v  €  Mj, 


J 
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where  ( v)j  denotes  the  inner  product  in  My  Let  be  the  Kronecker  delta.  Without  loss 
of  generality, 

0  <  <  ...  <  Xgj  S  Ay 

(Vf,4j))j  “  6ik’  1  <  '>k  <  Nj>  and  M? 

■jWPty?)  = 

With  each  space  My  we  define  discrete  norms 

IIMIlJ  =  fi-  ■=i0'!i,)‘,  for  v  =  £ 

i=  1  1=1 

where  we  have  suppressed  the  j  subscript  on  the  norm  and  — 2  <  s  <  2.  Note  that 
IIMIIj  —  IIMII  is  the  usual  energy  norm  on  level  j.  Hereafter,  we  drop  both  the 
superscripts  from  the  eigenvalues  and  eigenfunctions  (eigenvectors)  and  the  subscript  from 
the  dimension  of  the  spaces. 

We  require  the  bilinear  forms  on  adjacent  spaces  to  have  a  consistency  relationship. 
As  a  consequence,  the  energy  norms  on  two  adjacent  spaces  are  uniformly  consistent. 

Hypothesis  1.3  (Energy  Norm  Consistency):  For  all  j  >  1  there  exists  a  positive  constant 
Cj,  independent  of  j,  such  that 

aj(EjV,E.w)  =  Cj  a._j(v,w),  for  all  v,w  6  M}_y 

Since  |||v|||2  =*  aj(v,v).  This  hypothesis  implies  that  for  any  v  6  X._j,  the  energy  norm 
of  v  on  level  j— 1  is  equal  to  a  constant  times  the  energy  norm  of  E^v  on  level  j. 

We  now  define  and  analyze  a  k-Ievel  iterative  procedure  for  solving  Problem  1.2.  The 

process  involves  solving  problems  like  Problem  1.2  sequentially  for  j=l,  2 .  k.  The 

k-level  scheme  has  three  parameters:  m  and  n,  which  determine  the  number  of  smoothing 
iterations  used;  and  p,  which  is  used  in  a  recursion  iteration. 
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Algorithm  1.4:  MG(k,  m,  n,  p) 

Given  an  integer  k  >  0  and  {Mj,  a.(-,-),  fj(*)}]L,,  we  want  10  approximate  uk  6  Mk, 
where  ak(uk,v)  *»  fk(v)  for  all  v  6  Xk. 

(a)  If  k  =  1,  then  solve  directly. 

(b)  If  k  >  1,  then  one  iteration  of  the  k-level  scheme  takes  an  initial  guess  zQ  €  .Mk 
to  a  final  approximation  zm+n+i  €  Wk  in  three  steps: 

(i)  if  n  >  0,  define  Zj,  1  <  i  <  n,  by 

(*rVl>v)k  =  ^kHfk(v)  ~  *k(*i_i»v)l.  for  al1  v  ^  Mk.  (1.7) 

(ii)  Let  q  £  Mk-1  be  the  approximation  of  q  6  Mk_|  obtained  by  applying  p  iterations  of 
the  (k— l)-level  scheme  to  the  residual  equation 

ak_1(q,v)  -  C-1^^)  -  ak(VEkv)}  (1.8) 

=  fk-l(v),  for  all  v  €  -Mk_,, 

starting  from  an  initial  guess  zero.  Then  set 

Zn+1  "  Zn  +  Ek*  (L9) 

(iii)  If  m  >  0,  then  define  z.,  n+2  <  i  <  m+n+1,  by  (1.7). 

In  the  correction  recursion  iteration  (step  (ii)),  we  approximately  compute  the  elliptic 
projection  of  the  error  in  Mk_j  using  p  iterations  of  the  (k— l)-level  scheme  applied  to  a 
problem  of  Problem  1.2’s  form  with  j  *»  k  -  1.  In  the  smoothing  iterations  (steps  (i) 
and  (iii)),  error  components  whose  oscillation  are  “large”  are  damped.  A  simultaneous 
displacement  procedure  is  used  in  this  step.  Later  in  this  section  we  will  see  that  Ak  1  can 
be  replaced  by  a  particular  type  of  bound  (see  Hypothesis  1.0).  We  will  also  see  that 
(1.7)  can  be  replaced  by  other  iterations  which  are  computationally  more  attractive,  but 
do  not  affect  the  character  of  our  convergence  results.  The  use  of  iteration  (1.7)  simplifies 
the  initial  analysis  of  the  convergence.  Figure  1-2  contains  a  three-level,  two-iteration 
example  of  Algorithm  1.4  with  p  =*  2.  Note  that  computation  begins  with  the  largest 
space  and  uses  the  smaller  ones  only  to  solve  correction  recursion  problems. 


There  are  three  cases  of  note  in  Algorithm  1.4:  (a)  when  n  >  0  and  m  »*=  0,  (b)  when 
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Figure  1-2:  Three-Level  Example  of  Algorithm  1.4 

Two  iterations  on  level  three,  p  =  2 
Level 

1  ds  ds  ds  ds 

/  \  /  \  /  \  /  \ 

2  n  m+n  m  n  m*n  m 

t  \  /  Ai 

3  n  m-*n  m 

ds  =  direct  solve 

n.m  =  number  of  smoothing  iterations 


d=»0  and  m  >  0,  and  (c)  when  n  >  0  and  m  >  0.  Case  (a)  is  the  scheme  analyzed  by 
Astrakhantsev  [3],  Bank  and  Dupont  [6],  Hackbusch  [18,  19,  20],  Nicolaides  [24,  25],  and 
Van  Rosendale  [29]  for  finite-element  discretizations  of  various  elliptic  boundary-value 
problems.  Federenko  [14,  15]  and  Bakhvalov  [5]  analyzed  this  case  for  finite-difference 
discretizations.  Each  of  these  authors  assumed  in  their  convergence  proofs  that  p  >  1. 
Our  convergence  result  is  true  even  when  p  **  1.  Brandt  [10]  gave  a  heuristic  analysis  of 
(b)  and  (c)  for  finite-difference  discretizations  using  local  mode  analysis.  The  motivation 
for  studying  cases  (b)  and  (c)  comes  from  trying  to  understand  the  behavior  of  a  large 
finite-difference  program  which  is  described  in  Chapter  6  of  Douglas  [12,  13].  It  was 
observed  empirically  that  case  (b)  sometimes  required  fewer  correction  recursions  to 
achieve  the  same  accuracy  as  case  (a). 

Before  proving  a  convergence  theorem  for  Algorithm  1.4,  we  need  to  state  one 
definition  and  two  more  hypotheses.  The  first  hypothesis  is  a  bound  for  the  largest 
eigenvalues  and  the  other  is  an  error  estimate  for  the  correction  produced  by  the 
(k— l)-level  iteration.  Finally,  we  prove  a  lemma  describing  the  smoothing  iteration’s 
effect  on  the  error. 

Definition  1.5:  The  error  on  level  k  at  the  ith  stage  of  Algorithm  1 .4  will  be  denoted  by 


The  first  hypothesis  states  what  form  the  bound  for  the  maximum  eigenvalue  is 


assumed  to  hare. 
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Hypothesis  1.6  (Maximum  Eigenvalue):  There  exist  positive  constants  6  and  C,,,  each 
independent  of  j,  such  that 

Aj  <  C2N?5,  1  <  j  <  k. 

The  use  of  A k  in  the  smoothing  iteration  (1.7)  may  be  replaced  by  any  upper  bound 
satisfying  Hypothesis  1.6.  The  last  hypothesis  is  a  norm  estimate: 

Hypothesis  1.7  (Approximating  Error  Estimate):  For  some  a  with  0  <  a  <  1,  there 
exists  a  positive  constant  C3  such  that 

IIIEk<-'„IIU  <  CJNk2‘“lll«„lll1+„. 

For  problems  derived  from  elliptic  boundary-value  problems,  the  value  of  a  depends  on 
the  spatial  domain. 

The  following  lemma  is  used  in  the  convergence  proof  to  analyze  the  effect  of  the 
smoothing  on  elements  in  the  solution  spaces. 

Lemma  1.8:  Let  n  >  0  be  any  integer  and  zQ  6  Then  the  smoothing  iteration 

(1.7)  is  a  contraction  operator: 

illejll  <  IIMI-  (U°) 

Further,  for  every  fixed  0  <  oj  <  2  and  0  <  a  <  1, 

IIMI.  <  C5%J‘(2»  +  IlleJIU,-  (I  ”) 

Proof:  Let  N  ass  and  A  ss  yl^.  From  (1.7)  we  can  deduce  that 

(«-«i_,,Wk  —  ~  A~\(e.x_vr),  for  all  v  €  Affc,  i  >  1.  (1.12) 


We  can  expand  eQ  in  terms  of  the  eigenfunctions: 
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:n 

h 

BQ 


Using  (1.12)  we  can  show  that 

en  -  £  0;  (1  -  X;M)n  V»J.  (1.13) 

i=x 

Since  (1  —  X./yl)n  <  1,  we  have  that 

llle.HI  <  ll|e0ll|. 


The  proof  of  (1.11)  uses  (1.10)  and  the  Maximum  Eigenvalue  Hypothesis: 

hmc  =  ^xr<'  -  V'"2" 

=  //fyr%UAf  (i  -  xjAf' 

i=l  1  1  1  1 

<  Aa  max  |x°(l  -  x)2n|  £ 

*eio,i|  i=i  1  * 


<  Aa  (2n  +  a)-a  |||e0|||2_a 

<  qN2^(2n  +  ara|||e0||i2_o. 

Taking  the  square  root  of  both  sides  of  this  inequality  completes  the  proof. 

QED 

The  convergence  of  Algorithm  1.4  is  established  in  the  Theorems  1.9  and  1.10.  In 
essence,  these  results  say  that  the  error  on  level  k  can  be  reduced  by  any  positive  constant 
less  than  one  provided  the  correction  recursion  problem  on  level  k— 1  can  be  solved 
sufficiently  accurately.  Theorem  1.9  requires  that  p  >  1  and  is  true  for  any  k.  Theorem 
1.10  requires  that  k  be  finite  and  is  true  for  any  p. 
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Theorem  1.9  (Convergence  of  Algorithm  1.4):  Assume  that  Hypotheses  1.3,  1.6,  and 
1.7  hold.  Let  p  >  1  be  any  fixed  integer.  For  any  constant  0  <  7  <  1  there  exist  a 
nonnegative  integer  I  which  depends  only  on  p  and  7,  such  that 

Hlw,lll  <  7|||e„l!l,  for  all  m+n  >  I.  (1.14) 

Proof:  This  proof  is  motivated  by  the  work  of  Bank  and  Dupont  (6].  The  basic  idea  of 
this  proof  is  to  show  that  the  smoothing  iteration  (1.7)  reduces  the  oscillatory  components 
of  the  error  (corresponding  to  the  larger  eigenvalues)  while  the  correction  recursion 
iteration  (1.8)  reduces  the  smoother  components  of  the  error.  The  proof  is  by  induction 
on  the  index  k  of  the  space.  Assume  the  result  is  true  for  1,  2,  ...  ,  k— 1.  We  now  prove 
the  result  for 

By  Lemma  1.8  (with  w  =  1  +  a), 
llkjll  <  |||e0||| 

and 

ll|en|||1+0  <  C«/2Nf(2n  +  a)-<I/2|||e0|||.  (1.15) 

We  can  estimate  the  effect  of  the  correction  recursion  iteration  (1.8).  Using  (1.8)  and 
the  Energy  Norm  Consistency  Hypothesis,  we  have  for  all  v  6  Mk_j, 

ak_i(q-v)  —  crlak(vM 

or 

MEkcl_en-Ekv)  “  °-  (116) 

This  shows  that  the  exact  solution  q  of  the  (k— l)-!evel  correction  recursion  problem 
corrects  exactly  all  components  in  zn  which  belong  to  Take  v  =  q  in  (1.16).  Then 
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ak(Ekq-en,Ekq)  —  0  m  a^q^q)  -  ak(en,Ekq) 

=>  IIMtl  <  lllcJH  <  |||c0|||. 

Using  this,  the  Energy  Norm  Consistency  Hypothesis,  and  the  induction  hypothesis  gives 
us 


ll|Ek(q-q)|||  -  C71/2„,q_q,||  <  C^Vfllqfll  -  7p|||Ekq||| 
<  7p||ie0|||. 


(1.17) 


We  are  now  ready  to  estimate  |||em+n+1|||.  We  can  define 
en+l  "  E^-e„ 

=  (Ek^“en)  + 

If  Sm  reflects  the  effect  of  m  smoothing  iterations  on  any  v  g  .Mk,  then 

p  — 

cm+n+l  °  Cn+T 

Define  C4  s  C^C"  Using  Lemma  1.8  (with  u  =*=  1),  (1.15),  (1.17),  and  the 
Approximating  Error  Estimate  Hypothesis  yields 

lllem+n+l!H  <  ll|Sm(Ekq-en)|||  +  ||!SmEk(q-i)|!| 

<  C-/2(2m+a)-/2N^,|(Ekq_en„,i_a  +  j||Ek(q-q)||j 

<  q/2C»(2m  +  QrQ/2Nk^||!enjiiI+a  +  7p|||e0||| 

<  [C4(2m  +  ay°l\ 2n  +  a)~ °/2  +  7P)  |||e0|||.  (1.18) 

Choose  I  such  that 


C4(2m  +  a)-Q,/2(2n  +  a)~a!~  <  7  —  7P,  for  all  m+n  >  I. 


(1.19) 


Then 
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C 

I  ’ 

« 

;« 

X 


1 


I 


I 


« 


cm+n+l* 


S  7 


Clearly,  7  and  I  can  be  chosen  independent  of  the  N^. 


(1.20) 


QED 


When  p  =»  1,  the  proof  of  Theorem  1.9  breaks  down  because  the  right  hand  side  of 
(1.19)  is  zero.  If  we  assume  that  k  is  a  small  natural  number,  then  we  can  define  a 
sequence  {7j}j^j  such  that 

7j  -  7?.,  -  C5  >  0,  2  <  j  <  k  and  p  >  0, 

7j  >  0,  and 
7k  =  7<  1. 

Corresponding  to  Theorem  1 .9  is 


Theorem  1.10:  Assume  that  Hypotheses  1.3,  1.6,  and  1.7.  Let  p  >  0  be  any  fixed 
integer.  For  any  constant  0  <  7  <  1  there  exist  a  nonnegative  integer  I  which  depends 
only  on  p,  7,  and  C5,  such  that 

lllem+n+1HI  ^  7|||e0|i|,  for  all  m+n  >  I. 


Proof:  The  proof  follows  the  one  for  Theorem  1.9.  Equation  (1.17)  becomes 
ll|Ek(q-q)|||  <  7k_1|||e0|||. 
and  (1.19)  becomes 

C4(2m  +  Q)~0,/,2(2n  +  a)-0/2  <  7k  -  for  all  m+n  >  I. 


QED 

At  first  glance,  it  appears  that  choosing  m  =*  n  in  (1.18)  would  be  better  than 
choosing  either  n  =  0.  m  >  0  or  n  >  0,  m  =■  0.  However,  special  case  proofs  show  that 
this  is  not  necessarily  true.  Assume  that  n  >  0,  m  =  0.  By  modifying  the  proof  of 


Theorem  1.9,  we  can  show  that 
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IIKh.,111  <  (C>/2(2n  +  «)-/’  +  V!  |||«0|||. 

For  n  large,  this  bound  for  the  error  reduction  looks  like 
Cj/2  (2n)“a/2  <  €  m  7  -  -f. 


So 

n  >  lCi/0r^_2/o, 

—  2  4 

is  required  to  reduce  the  error  by  a  factor  of  t.  If  S'  =*  m,  the  error  reduction  (see  (1.18)) 
looks  like 

C4(2ST*  <  e. 

Once  again, 

S'  >  },C\'arl/°  mm  tx!a  -\CxJa  C*la  »  tx!a n 

is  required  to  reduce  the  error  by  a  factor  of  e.  The  amount  of  work,  which  depends  on  n 
and  2n^  depends  on  the  relative  sizes  of  one  and  2°€.  Thus,  we  have  shown 

Theorem  1.11:  Assume  Theorem  1.9  holds.  Define  t  =  7  —  7P.  Then  choosing  n  >  0, 
m  =  0  in  (1.18)  is  better  than  choosing  S'  =  m  >  0  only  when  2°t  >  1.  Alternately, 
choosing  5T=  m  >  0  in  (1.18)  is  better  only  when  2°<  <  1. 


The  practical  significance  of  this  is  that  if  7  —  7P  (or  -  qj|_|)  is  large,  it  is  more 
efficient  to  do  the  smoothing  at  once,  rather  than  splitting  it  around  the  correction 
recursion.  Similar  analysis  holds  for  the  case  of  n  =  0,  m  >  0. 

We  now  analyze  the  cost  of  one  iteration  of  Algorithm  1.4.  Let  F(N)  be  the  cost  of 
reducing  the  error  by  a  factor  of  7  for  a  problem  with  N  unknowns.  We  assume  that  the 
cost  of  the  smoothing  iterations  (1.7)  (or  an  iteration  with  similar  properties)  on  the  finest 
level,  level  k,  can  be  bounded  by  C6(m+n)Nk  ss  C7Nk,  where  Cg  is  independent  of  k. 
The  cost  of  the  correction  recursion  (1.8)  is  pF(Nk_j).  Thus, 


F(Nk)  -  PF(Nk_,)  +  C7 Nk. 


(1.21) 


PC 


5  ' 


< 


N 


ft 


4 


( 


Since  Nk  ~  <rNk_j  (see  (1.2)),  the  solution  of  (1.21)  is 

F(Nk)  -  C7Nk  h  (p/4- 

Asymptotically  in  k, 

CftN,  1  <  P  <  o 


F(N)  <  ^  CgNlogN,  p  =*  a 
CftN,0SP,  p  >  ctj 


(1.22) 


where  the  logarithms  are  base  a  [2].  One  of  the  surprises  of  multi-grid  is  that 
implementations  (e.g.,  Bank  and  Sherman  [7]  or  Douglas  [12,  13])  exhibit  the  asymptotic 
convergence  rates  for  very  small  k’s. 

Choosing  1  <  p  <  o  leads  to  an  optimal  order  algorithm,  in  the  sense  that  the  error 
can  be  reduce  by  a  fixed  factor  of  7  each  iteration  with  work  proportional  to  the  number 
of  unknowns.  When  k  is  finite,  the  cost  for  the  direct  solution  of  the  one  level  problem 
may  not  be  majorized  by  C7Nj}.  In  this  case,  the  algorithm  is  not  optimal  order. 

We  may  want  to  reduce  the  initial  error  by  a  factor  of  N-q^or  some  fixed  q.  The 
obvious  implementation  would  then  require  F(N)logN  operations.  We  assume  the 
solutions  u-  of  Problem  1.2  satisfy 

IllPjU— UjlH  <  KNpfj  >  1, 


where  K  is  a  constant  independent  of  Nj.  Denote  by  uj  the  computed  solution  of  the  j- 
level  scheme  (Algorithm  1.4).  To  avoid  the  extra  logN  factor,  we  use  Ejuj_j  as  the  initial 
guess  to  uj,  j  >  1,  and  prove  that  the  initial  error  is  small.  Approximate  solutions  uj  of 
finite-dimensional  Problem  1.2  are  generated  using  the  following: 


Algorithm  M2:  MGN(j,  r,  m,  n,  p) 

Given  an  integer  j  >  0  and  {M[t  a-(-,-),  f(  )}(—l,  we  want  to  approximate  u.  6  Wj,  where 
a.(u.,v)  =  f.(v)  for  all  v  €  .M:- 


Figure  1-3:  Three- Level  Example  of  Algorithm  1.12 


One  iteration  (r  =  1)  on  level  three,  p  =  2 
Level 

1  ds  ds  ds  ds 

w  \  /  \  /  \ 

2  n  m  n  m+n  m 

\  /  * 

3  n  m 

ds  =  direct  solve 

n,in  =  number  of  smoothing  iterations 


(a)  If  j  ==  1,  then  solve  Problem  1.2  directly. 

(b)  If  j  >  1,  then  starting  from  an  initial  guess  *0  as  E.uj_j,  apply  r  iterations  of 
MG(j,  m,  n,  p)  (see  Algorithm  1.4)  to  Problem  1.2  to  obtain  uj. 


This  algorithm  actually  has  four  parameters.  The  parameter  r  determines  the  number  of 
iterations  of  MG(j,m,n.p)  to  use.  Recall  that  m  and  n  are  the  number  of  smoothing 
iterations  and  that  p  is  the  number  of  correction  recursion  iterations  in  Algorithm  1.4. 

Figure  1-3  contains  a  three-level  example  of  one  iteration  (r  —  1)  of  Algorithm  1.12. 
For  the  correction  recursion  problems,  p  “  2.  Unlike  Algorithm  1.4,  computation  begins 
with  the  smallest  space  and  winds  its  way  “down"  to  the  largest  space.  It  is  worth  noting 
that  when  n  =  0  and  *=  0  in  Figure  1-2,  ,'ie  two  algorithms  become  much  more 
similar.  In  this  case,  the  First  recursion  correction  problem  (see  (1.8))  in  Figure  1-2  has 
f.  (v)  »  f.  (v). 


The  convergence  properties  of  Algorithm  1.12  are  stated  in  the  Theorem  1.13. 
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Theorem  1.13:  Assume  that  Hypotheses  1.3,  1.5,  and  1.7  hold.  Let  r  >  1  be  any  fixed 
integer.  Suppose 

(i)  !l|Pju— ujlll  <  KNpf  j  >  1, 

(«)  IIIPju— EjPj_1u|||  <  KNpf  j  >  2,  and 

(m)  lll»rSi|||  <  KNp/ 

Then 

(a)  For  any  integer  p  >  1,  there  exists  a  nonnegative  integer  I  such  that  for  j  >  1, 

IllUj— uj|l|  <  KNpf  for  all  m+n  >  I. 

(b)  HlPjU-ujlll  <  2KNpf 

(c)  the  cost  of  computing  uj  is  bounded  by  CfiF(Nj),  where 
Cfl  as  <rr/(«r— 1)  is  independent  of  j. 

Proof:  (a)  The  proof  is  by  induction  on  the  index  j  of  the  space.  Assume  the  result  is  true 
for  1,  2,  ...  ,  j— 1.  We  now  prove  the  result  for  Mj.  Define  ej  s  |||uj — 3j|||  for  j  >  1. 
After  r  iterations  of  the  j-level  scheme, 

ej  <  V  !l|urEj3j_1||| 

<  "f  {  lll«j-Pjulll  +  IIIPj«-Ejpj_,u|H  + 

ll|Ej(pj_iU-uj_1)|||  +  HlEj(uj_i— uj_l)U|  } 

<  V  {  K(2  +  C1/2  +  Cp^)Npi 

Choose  7  <  (2  +  C|/2  +  Cp/3,q^)~1/r  and  I  according  to  Theorem  1.9  or  Theorem  1.10. 
Hence, 
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(b)  This  is  a  simple  consequence  of  the  triangle  inequality: 

IllPjU-Srill  <  IIIpjU— Ujlll  +  |||uruj||| 

<  2KNp* 

(c)  The  cost  of  computing  uj  is  bounded  by 

F(Nj)  +  r  F(Nj)  <  Cg  F(N-)  or/(<r-l)  =  Cg  F(Nj). 

QED 

It  is  worth  pointing  out  that  r  is  independent  of  j.  This  tells  us  that  Algorithm  1.12  is 
optimal  order  whenever  1  <  p  <  cr. 

In  conclusion,  we  have  shown  that  solutions  to  problems  like  Problem  1.1  can  be 
approximated  in  finite-dimensional  spaces  using  an  optimal  order  procedure.  In  the  next 
section  we  verify  that  Hypotheses  1.3,  1.6,  and  1.7  hold  for  particular  discretizations  of 
several  elliptic  boundary-value  problems. 


2:  Applications 


2.1  Model  Problems 

In  this  section,  we  discuss  the  application  of  the  algorithms  of  Section  1  to  the 
solution  of  large  spane  linear  systems  which  arise  in  connection  with  finite-element  and 
finite-difference  procedures  for  solving  self-adjoint  elliptic  boundary-value  problems.  We 
show  a  few  examples  in  which  the  algorithms  and  theorems  of  Section  1  apply.  We  will 
see  that  these  can  be  optimal  order  results. 

Section  2.2  is  concerned  with  a  finite-element  procedure.  Our  model  problem  is  the 
Neumann  problem 


—  v(PV«)  +  Su  =  f  in  f} 


u  =■  0  on  dfl, 
n 


where  Cl  is  a  polygonal  domain  in  R2.  We  assume  that  P  £  C^/J),  S  €  and  that 

there  exist  positive  constants  p,  p,  s,  and  s  such  that 


0  <  p  <  P(x)  <  p  and  0  <  s  <  S(x)  <  s,  for  all  x£/7. 


Most  of  our  arguments  apply  to  the  Dirichlet  problem 


—  v(PVu)  +  Su  =  f  in  Cl 


u  =>  0  on  3j? 


with  only  minor  modifications.  We  comment  on  the  extensions  as  they  arise. 

In  Section  2.3,  we  discuss  finite-difference  approximations  to  the  model  problems 


-  (Pux)x  +  Su  —  f  in/?  m  (0,1) 
u(0)  —  u(I)  —  0 


and 
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"  (Pux)x  -  (PuA  +  Su  "  fin«-(0,lW0,l) 


u  as  0  on  dfl. 


We  assume  that  f  €  L2(/7)  and  that  P  and  S  are  constrained  as  before. 


2.2  Finite-Elements 

We  seek  a  -weak  form  solution  of  (2.1):  find  u  £  1/*(J?)  «uch  that 

a(u,v)  =>  (f,v)  for  all  v  £  #*(/?),  (2.6) 

where 

a(u,v)  =  /  Pyu-Vv  +  Suv  dx  and 
n 

(2.7) 

(f,v)  ■  /fvdx. 

J7 

Then  there  exists  a  unique  weak  solution  u  £  )/*(/?)  for  all  f  £  L2(l?)  (see  Ciarlet  (l  1]). 
The  spaces  l/s,  for  s  a  positive  integer,  will  be  the  usual  Sobolev  spaces  equipped  with 
norms 

Bull?  =  S  IJD^uUj  -  E  (D^u,D^u). 
s  w<»  w<* 

The  spaces  H3  for  s  positive  and  non-integral  will  be  defined  by  interpolation  (see 
Agmon  [l]  or  Lions  and  Magenes  [23]).  For  s  negative,  will  be  defined  as  the  dual  of 
l/-s.  The  bilinear  form  a(-,  )  induces  the  energy  norm  |||u|i|2  =  a(u,u). 

A  modest  amount  of  elliptic  regularity  for  the  solution  of  (2.6)  is  required. 

Hypothesis  2.1  (Regularity):  We  assume  there  exists  a  constant  0  <  or  <  1,  such  that  for 
all  f  €  there  exists  a  unique  solution  u  £  )t1+0  of  (2.6)  and 

Ml+a  <  c(p ,s,n)  ||f||a_I. 
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For  a  complete  discussion  of  what  values  of  a  correspond  to  specific  domains  1 1 ,  see 


Kellog  [22]  and  Babuska  and  Aziz  [4]. 

We  now  consider  a  finite-element  approximation  of  (2.6).  Let  Tj,  j  >  1,  be  a  nested 
sequence  of  triangulations  of  fj.  Take  Tj  to  be  a  fixed  triangulation.  For  T  6  Tj,  denote 
the  diameter  of  T  by  hp,  and  let  iLp-dp  denote  the  diameter  of  the  inscribing  circle  for  T. 
Define 


h.  as  maxh-v 

1  T6T,  1 


and 


The  constant  $0  is  a  measure  of  the  shape  regularity  of  the  triangles  in  Tj  and  is  a 
measure  of  their  uniformity.  We  construct  Tj,  j  >  1,  inductively:  divide  every  T  €  T-_j 
into  congruent  triangles,  where  is  independent  of  j.  When  p  =  2,  this  means  we 
construct  four  triangles  in  Tj  by  pairwise  connecting  the  midpoints  of  the  edges.  Each 
triangulation  Tj  will  have  shape  regularity  and  uniformity  constants  SQ  and  and  will 

have  h-  ss  maxh-  =*  yi1-Jh.. 

J  TeTj  1  1 

With  each  triangulation  Tj,  we  associate  the  Nj-dimensional  space  M-  of  C°-piecewise 
linear  polynomials.  Following  (1.2),  we  know  that 

Njtl  ~  «Nj.  (2.8) 

where  <t  =  n“  asymptotically.  Since  the  triangulations  are  nested,  we  have  that  Mj  is  a 
subset  of  Mj+1,  j  >  1.  The  spaces  Mj  satisfy  the  following  standard  approximation 
property  [8,  9,  21,  27]:  if  u  €  1/*,  1  <  s  <  1+a,  then  there  exists  a  Uj  €  Mj  such  that 

||u-Uj|]0  +  hjllu-Ujll,  <  C( S0,6vn)  h?  Bull,.  (2.9) 

We  briefly  remark  on  the  Dirichlet  model  problem  (2.3).  The  definition  of  Tj,  j  >  1, 
remains  the  same.  Let  Mj  be  a  subset  of  Mq  and  let  it  be  the  space  of  C°-piecewise  linear 
polynomial  associated  with  Tj  satisfying  the  Dirichlet  boundary  conditions.  Then  Mj  is  a 
subset  of  Mj+1,  j  >  1,  as  before.  With  this  modification,  the  results  of  this  section  will 
remain  valid. 

Using  the  notation  of  Section  1,  we  define  for  each  space  Mj,  j  >  1, 
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aj(u,v)  m  a(u,v),  f^(v)  m  (f,v),  and  (u,v)j  as  (u,v)  for  all  u,v  €  My 

The  interpolation  and  projection  operators,  Ej,  Pj,  and  ij  (see  (1.3)  and  (1.4)),  are  the 
natural  projections  and  injections. 

Associated  with  M-  are  eigenvalues  \p),  eigenfunctions  y>p\  and  a  maximum 
eigenvalue  .4-  of  a(v)  fulfilling  the  requirements  of  (1.6),  where  1  <  i  <  N-.  A  simple 
homogeneity  argument  shows  that 

dj  <  C(P,S ,s0,6vn)h-2  (2.10) 

(see  Strang  and  Fix  (28j).  For  — 2  <  s  <  2,  we  define  discrete  norms 


(2.11) 


where  we  have  suppressed  the  j  subscript  on  the  norm.  Note  that  |||v|||j  **  |]|v|jj  and 
|||v|||0  is  comparable  to  ||v||0.  In  fact,  the  proof  of  the  following  norm  equivalence  is 
almost  identical  to  Lemma  1  in  Bank  and  Dupont  [6]: 

Lemma  2.2:  There  exists  a  constant  C  =  CtPjS,/?,^,^,/?)  such  that  for  0  <  s  <  1, 
C-l||v||s  <  |||v|!!s  <  C||v||s. 


In  order  to  establish  the  convergence  of  Algorithms  1.4  and  1.12  for  the  finite*element 
case,  we  must  verify  Hypotheses  1.3,  1.6,  and  1.7.  It  is  immediate  that  (the  Energy  Norm 
Consistency)  Hypothesis  1.3  holds.  Using  (2.10)  and  the  fact  that  Nj  ~  Chr2we  can 
verify  that  (the  Maximum  Eigenvalue)  Hypothesis  1.6  holds. 

A  duality  argument  is  used  to  verify  (the  Approximating  Error  Estimate)  Hypothesis 
1.7,  When  a  is  an  integer,  this  is  a  standard  result  (11,  28}. 

Lemma  2.3:  Let  a  be  defined  by  the  Regularity  Hypothesis.  Then 


lllq-*nIIU  <  C,2ohk2o,|||en!||1+o. 


Proof:  By  Lemma  2.2  and  duality, 
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llfi-.HI,-.  <  cii-J,.. 

(p,q-en) 

=»  C  sup - 

pe^-1 

For  p  €  let  rj  6  ^*+l  be  defined  by 
a(»j,v)  =—  (p,v)  for  all  v  €  #l. 

Taking  v  —  q-en  gives  us 
(P.q-en)  —  a(«j,q-en) 

=  a(»j— w,q— en),  for  any  u  6 
By  the  Regularity  Hypothesis  and  (2.9), 

(m-,)  <  a*.,  mu,  nfi-.ni 

<  <V*  hf  114.-,  Ilfi— .III- 

Combining  this  with  (2.12)  yields 

l!lq-«nlll,^  <  Cm*  h£  ||lq-en|||. 

However, 

lllq-«nlll2  -  afq-en.q-en) 

“  -  a(i-en,en) 

Substituting  (2.13)  into  the  right-hand  side  of  (2.14)  gives  us 
lllq-enll!  <  CM*h*|||en|||I+0, 


(2.12) 


(2.13) 


(2.14) 
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which  we  substitute  back  into  (2.13)  to  complete  the  proof. 

QED 

This  proves  that  Algorithm  1.4  converges  at  the  rates  specified  by  Theorems  1.9  and 
1.10  for  the  finite-element  case  of  this  section. 

One  of  the  advantages  of  finite-element  methods  is  that  the  theory  of  Section  1  can 
be  applied  using  a  variety  of  norms.  As  an  example,  we  prove  a  special  case  of 
Theorem  1.9  for  the  L2  norm.  It  is  similar  to  the  results  of  Nicolaides  [25]  for  the  l2  norm 
and  is  the  analogue  of  Corollary  1  of  Bank  and  Dupont  [6].  To  get  L2  results,  we  assume 
that  the  solution  u  has  ^2  regularity,  i.e.,  a  =  1.  This  assumption  requires  that  Cl  be 
convex  [17,  28]. 

Theorem  2.4:  Assume  the  Regularity  Hypothesis  holds  for  a  =  1.  Let  p  >  1  be  any 
fixed  integer.  For  any  constant  0  <  7  <  1  there  exists  a  nonnegative  integer  I,  which 
depends  only  on  p  and  7,  such  that 

Hem+n+lflo  ^  T'IMo  for  3,1  m+n  ^  1 

Proof:  The  proof  is  by  induction  on  the  index  k  of  the  space.  Assume  the  result  is  true 
for  1.  2,  ...  ,  k— 1.  We  now  prove  the  result  for  .M^.  From  (1.13)  it  is  immediate  that 

kilo  <  kilo-  (2-15) 

Using  an  argument  similar  to  the  proof  of  Lemma  1.8  shows  that 

kll2  -  C^-2hk?‘2n+ir1||e0||0.  (2.16) 

Lemmas  2.2  and  2.3  with  a  =  1  yield 

lh-«„llo  <  Cx*  h;  lie, J|2.  (2.17) 

The  analogue  of  (1.17)  is  derived  using  the  induction  hypothesis  and  (2.15)  —  (2.17): 
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llq-illn  <  ^PMIC 


<  ^{liq-eJo  +  lleJIo) 


<  7p{C/i2h2||en||2  +  |je0|j0} 


<  7p{C(2n+l)-1  +  l}||efl|L 


<  C7p||e0||0. 


Using  an  argument  similar  to  (1.18)  gives  us 


l|em+n+1ll0  <  C{(2m  +  l)-1(2m-l)-1  -»-  ||e0||0 


and  Theorem  2.4  follows. 


When  p  =  1  we  can  prove  the  analogue  to  Theorem  1.10. 

It  is  useful  to  associate  with  Xj  a  symmetric,  positive  definite  bilinear  form  b-(v), 
which  is  comparable  to  the  L2  inner  product.  We  assume  there  exists  a  constant  /?, 


independent  of  hj,  such  that 


0  <  d  V.v)  <  bj(v,v)  <  (i(v,v)  for  all  v  £  X.,  v  0. 


(2.18) 


For  bj(v),  we  define  generalized  eigenvalues  Xt  and  eigenfunctions  1  <  i  <  N^,  by 
a(v,Vj)  =  ^  b.(v,Vi)  for  all  v?^,  0  <  X'j  <  T2  <  ...  <  FN  sb  X, 

bj(v;.vk)  =  *ik,  and  a(?j,vk)  =*  1  <  k  <  N-. 


It  then  follows  that 


A  =  max 


a(v .  v) 


<  d  max 


a (v. v) 


v€j trCv.v)  v6Xj.»^0  (v.v) 


<  CCP.S.a^.^j.^h" 


i 
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For  — 2  <  s  <  2,  define  |||v|||2b  similarly  to  |||v|||2  (see  (2.11)).  Note  that 
|| jv||J j  b  =  || jv|[|  and  |j|v|||ob  is  comparable  to  ||v||0.  In  fact,  a  norm  equivalence  similar 
to  that  in  Lemma  2.2  can  be  proved  ■with  |i[v|||gb  substituted  for  |||v|||s. 

The  smoothing  iteration  (1.7), 

(*-Vl’v)k  =  ^ 1  (fk(v)  -  \(^vy)],  ^  all  v  6  Mk, 

requires  the  solution  of  a  linear  system  involving  the  mass  matrix  at  every  step.  In 
practice  this  is  too  expensive.  We  can  replace  (1.7)  with  the  smoothing  iteration 

bk(zi_zi-l’v)  =  Kf>v)  ~  a(z;-i'v)l»  for  a11  v  ^  -Mk-  (2.19) 

There  are  numerous  choices  for  b-(v),  j  >  1.  When  the  standard  nodal  basis  is  used, 
(2.18)  is  satisfied  by  b(v)  corresponding  to  the  diagonal  of  the  mass  matrix.  This  allows 
smoothing  by  an  under-relaxed  Jacobi  scheme. 

The  convergence  of  the  k-level  scheme  (Algorithm  1.4)  is  summarized  by  the 
following  result.  Its  proof  is  almost  identical  to  the  ones  for  Theorems  1.9  and  2.4. 

Theorem  2.5:  Assume  the  Regularity  Hypothesis  holds.  Define  the  k-level  scheme 
(Algorithm  1.4)  using  (2.19)  instead  of  (1.7).  Let  ||-||  denote  either  the  energy  or  L2  norm 
(and  be  fixed).  Let  p  >  1  be  any  fixed  integer.  For  any  constant  0  <  y  <  1  there  exists 
a  nonnegative  integer  I,  which  depends  only  on  p  and  y,  such  that 

llem+n+lll  <  7l|e0||  for  all  m+n  >  I. 

We  can  extend  this  theorem  to  cover  the  case  of  p  =  1  analagously  to  Theorem  1.10. 

We  conclude  this  section  by  noting  that  for  the  finite-element  method  of  this  section, 
Algorithm  1.12  converges  at  the  rate  specified  by  Theorem  1.13  for  either  the  energy  or 
L2  norm.  Let  ||  ||  denote  either  of  these  norms.  Then  Theorem  1.13  can  be  rewritten  as 


follows: 
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Theorem  2.6:  Assume  the  Regularity  Hypothesis  holds.  Let  r  >  1  be  any  fixed  integer. 
Suppose 

(i)  ||u— UjH  <  KhJ,j  >  1,  and 

(ii)  Huj-uJI  <  Kh}. 

Then  for  any  integer  p  >  1,  there  exists  a  nonnegative  integer  I  such  that  for  j  >  1, 

llu- — ut|l  <  Kh?for  all  m+n  >  I. 

"  j  j"  —  J  _ 

Moreover, 

l|u-3j||  <  2Khj>, 

and  the  cost  of  computing  uj  is  bounded  by  CflF(Nj),  where 
Cfl  =  crr/(<r —  1 )  is  independent  of  j  and  a  is  defined  by  (2.8). 

It  is  worth  pointing  out  that  r  is  independent  of  j.  This  implies  that  Algorithm  1.12  is 
optimal  order  when  1  <  p  <  or- 

2.3  Finite-Differences 

As  in  Section  2.2,  we  seek  weak  form  solutions  to  both  (2.4)  and  (2.5).  We  prove  that 
the  theorems  and  analysis  of  Section  1  apply  to  a  particular  discretization  of  (2.4)  and 
(2.5).  Most  of  this  section  consists  of  an  analysis  of  a  constant  coefficient  problem.  In 
Section  2.3.3,  we  explain  how  to  extend  the  analysis  to  the  variable  coefficient  self-adjoint 
problems  (2.4)  and  (2.5). 

2.3.1  One-Dimensional  Definitions 

Let  k  >  1  and  NQ  >  0  be  fixed  integers.  We  define  Nj  s  2N-_j  +  1,  j  >  1. 
Following  (1.2),  we  know  that  <7  =  2.  We  define  k  uniform  grids  /?*,  1  <  j  <  k,  by 

hk  =  (Nk  +  I)"*,  h.  s  2k~\,  and  a  {  ih.  |  1  <  i  <  N.  }. 


For  the  moment,  we  assume  that  Psl  and  S  s  0  in  (2.4).  The  general  case  will  be 


analyzed  at  the  end  of  the  section.  We  discretize  (2.4)  by  approximating  the  second 
derivatives  by  central  differences  to  obtain  a  system  of  linear  equations 


Akuk  —  hkfk’ 


(2.20) 


where  uk  and  fk  approximate  the  solution  u  and  the  right-hand  side  f  on  rf*.  The 
Nk  x  Nk  matrix  Ak  is  given  by 


2  -1 

-1  2  -1 

-1  2 


(2.21) 


■1  2  -1 
-1  2 
-1 


(see  Varga  [30]).  Factoring  and  solving  this  algebraic  system  directly  requires  only  5Nk 
multiplications  and  3Nk  additions  [16].  We  will  show  that  the  multi-grid  algorithms  of 
Section  1  require  0(NklogNk)  operations  to  solve  this  problem  to  accuracy  comparable  to 
the  discretization  error,  i.e.,  0(Nk^. 

We  associate  a  solution  space  At-,  1  <  j  <  k,  with  each  grid.  Let  N  =  Nk.  Then 
Mk  =5  (  v  |  v  6  RN  }  and 

A^k _ j  53  {  v  I  V  €  Mk  and  (2.22) 

v  =  <  0,  vr  0,  v2,  0,  ...  ,  0,  vN^  ,  0  >  }. 

The  spaces  .Mk_2,  ...  ,  are  defined  recursively. 

We  define  picecwise-linear  interpolation  operators  between  adjacent  solution  spaces 


E2:  Mj-i  “*  and  R: 


by  means  of  the  matrices 
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0  .5  0  0 

0  1 

0  .5  0  .5 

1 

0  0  .5 


E 


2 


0 


0  0 


0  0 
0  .5 


0  0 


(2.23) 


1 

0  0  .5  0  .5  0 

1  0 
0  0  .5  0 


and  R  =  ET. 

Problems  on  coarser  grids  are  defined  using  the  interpolation  matrices.  Linear 
systems  on  coarser  grids,  Aj,  1  <  j  <  k,  are  defined  recursively: 


k— i 


E2  Ak-j+l  E2‘ 


Let  r  €  RN  be  a  right-hand  side  for  level  k.  Then  the  (k-J)-Uvel  problem  is  defined  by 


El  r  m 


fk-i-  *W 


(2.24) 


Half  the  rows  and  columns  of  Akl  are  zero.  We  can  re-order  the  matrix  so  that  the 
nonzero  rows  and  columns  are  the  first  Nk_j  rows  and  columns.  Then  Ak_j  has  a 
submatrix  whose  form  is  identical  to  Afc.  The  re-ordered  solution  space  has  the 

form 


{  v  6  ;Mk  |  v  —  <  Vj,  v2,  ...  ,  vNk  i  ,  0,  ...  ,  0  >  }. 

The  (k— i)-level  problems  are  defined  similarly. 

Before  the  theory  of  Section  1  can  be  applied,  we  must  complete  the  definitions 
required  for  the  triples  used  by  Problems  1.1  and  1.2.  For  j  >  1,  a  bilinear  form  aj(-,-),  a 
linear  functional  fj(v),  and  an  inner  product  (-,-)j  are  defined  by 

aj(u,v)  as  uT  Aj  v,  fj(v)  s  hj  fjv,  for  all  u,v  g  Mj,  and 
(u,v)j  =  ^  uT  v. 


(2.25) 
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The  bilinear  form  a.(v)  induces  the  energy  norm  |||u|||2  m  a-(u,u)  for  all  u  €  Mj. 
Associated  with  each  a^-,  )  are  N  a  Nj  nonzero  eigenvalues  x[^  and  eigenvectors 
satisfying  (1.6).  Let  Aj  be  the  largest  eigenvalue.  For  — 2  <  s  <  2,  discrete  norms  are 
defined  by 


(2.26) 


where  we  have  suppressed  the  j  subscript  on  the  norm.  Note  that  1 1 |v] 1 1 1  =  |||v|||  is  the 
usual  energy  norm  on  level  j.  Hereafter,  we  drop  the  superscripts  from  the  eigenvalues 
and  eigenvectors. 

2.3.2  Two-Dimensional  Definitions 

Let  k  >  1  and  NQ  >  0  be  fixed  integers.  We  define  Nj  ss  2N^_j  +  1,  j  >  1.  We 
define  k  uniform  grids  ,  1  <  j  <  k,  as  products  of  the  one-dimensional  domains  flh 

hk  3  (Nk  +  1)-*,  h.  m  2k-\,  and  m  ff  ®  fl». 

Each  grid  covers  the  interior  of  /?  with  Nj  as  N?  points.  Following  (1.2),  we  know  that 
<7  3  4. 

For  the  moment,  we  assume  that  P  3  1  and  S  3  0  in  (2.5).  The  general  case  will  be 
analyzed  in  Section  2.3.3.  We  discretize  (2.5)  by  central  differences  to  obtain  a  system  of 
linear  equations 

Bkuk  =  b*fk.  (2.27) 

where  uk  and  fk  approximate  the  solution  u  and  the  right-hand  side  f  on  ffy.  The  Nj  x  Nj 
matrix  Bk  is  given  by 

Bk  ~  \  (Ak  ®  *N  +  *N  ® 

where  Ij,  is  the  NkxNk  identity  matrix  and  Ak  is  defined  in  (2.21). 

We  define  solution  spaces  Mj,  1  <  j  <  k,  as  the  tensor- product  of  the  one¬ 
dimensional  solution  spaces  defined  in  (2.22).  Following  the  techniques  of  Section  2.3.1, 
we  define  piecewise-bilinear  interpolation  operators  between  adjacent  solution  spaces 


Problems  on  coarser  grids  are  defined  using  the  interpolation  matrices.  Linear 


systems  on  coarser  grids,  Bj,  1  <  j  <  k,  are  defined  recursively: 


B 


k-j  “  EbBk-j+lEb- 


Let  N  s  Nk  and  r  £  RN0RN  be  a  right-hand  side  for  level  k.  Then  the  (k—l)-level 
problem  is  defined  by 


Bk-i^  -  Ebr  =  fk-i>  Mk-r 


(2.28) 


Similar  to  the  one-dimensional  case,  we  define 

a.(u,v)  =  uT  B.  v,  f-(v)  =  h^fTv,  for  all  u,v  6  Mj,  and 

(u.v^  ss  h?uTv. 

The  bilinear  form  a.(-,  )  induces  the  energy  norm  |||u|||2  ss  a.(u,u)  for  all  u  6  .NL. 
Associated  with  each  aj(-,-)  are  Nj  nonzero  eigenvalues  \f/u  and  eigenvectors  i> ^  satisfying 
(1.6).  Let  A •  be  the  largest  eigenvalue.  For  —  2  <  s  <  s,  discrete  norms  |||v|||  are 
defined  similarly  to  (2.26). 


2.3.3  Convergence 

In  order  to  establish  the  convergence  of  Algorithms  1.4  and  1.12  for  the  finite- 
difference  cases  of  this  chapter,  we  must  verify  Hypotheses  1.3,  1.6,  and  1.7.  In  this 
section,  E  will  refer  to  either  E0  or  Eb.  We  begin  by  showing  that  (the  Energy  Norm 
Consistency)  Hypothesis  1.3  holds.  It  is  derived  using  the  definition  of  the  linear  systems 
A._j  and  Bj_1,  j  >  1. 

Lemma  2.7:  Let  j  >  1  be  an  integer.  Then 
a.(Ev,Ew)  =*  a._l(v,w),  for  all  v,w  £ 


That  (the  Maximum  Eigenvalue)  Hypothesis  1.6  holds  is  a  simple  consequence  of  the 


► 

f. 
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explicit  formula  for  the  eigenvalues  of  A-  and  B._j. 

Lemma  2.8:  /L  <  4dh~,2  where  d  is  the  dimension  of  the  problem. 

The  proof  that  Hypothesis  1.7  holds  is  a  simplified  version  of  the  proof  of  Lemma 
2.3  with  o  =  l.  It  requires  Lemma  2.9. 

Lemma  2.9:  Let  k  >  1  be  an  integer  and  d  the  dimension  of  the  probelm.  Let 
C(s,l)  =  4(s-1)/2  and  C(s,2)  =  8s/2.  For  s  =  0,  1,  or  2  and  u  6 

min  |||u-Ev|||s  <  C(s,d)  h2"s  |||n|||,.  (2.29) 

The  proof  is  contained  in  Douglas  [13].  This  proves  that  Algorithm  1.4  converges  at  the 
rates  specified  by  Theorems  1.9  and  1.10  for  the  finite-difference  cases  of  this  section.  As 
a  consequence  of  Lemma  2.9  we  can  estimate  the  constant  C4  in  (1.18),  which  governs  the 
number  of  smoothing  steps  required  in  Algorithm  1.4.  For  the  one-dimensional  problem 
(2.4),  C4  =  4  and  for  the  two-dimensional  problem  (2.4),  C4  =  16-21/2. 

Before  we  can  show  that  the  j-level  scheme  (Algorithm  1.12)  holds,  we  need  to  define 
projection  and  injection  operators  between  the  spaces  H  and  My  1  <  j  <  k  (see  (1.3)): 

Pj :  H  — — M-  as  evaluation  of  u  €  H  at  the  grid  points  of  CP  and 

'j:  H  as  piecewise-linear  interpolation. 

These  operators  have  properties  which  are  worth  pointing  out,  namely, 

p.i.  =  Identity  on  and  =  ij_jPj_i  on 

We  can  verify  that  the  assumptions  of  Theorem  1.13  hold  using  these  properties  and 
simple  known  facts  about,  the  model  problems.  Details  are  contained  in  Douglas  [13]. 
This  shows  that  Algorithm  1.12  converges  at  the  rate  specified  in  Theorem  1.13. 

We  conclude  this  section  by  considering  (2.4)  and  (2.5)  when  they  are  no  longer 
restricted  to  P  as  1  and  S  ss  0.  We  discretiie  (2.4)  and  (2.5)  by  central  differences  to  get 
Nk  x  Nk  linear  system  of  equations.  As  before,  we  define  coarse  grid  approximations  of 
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the  problems  using  the  interpolation  matrices.  We  define  inner  products  ST(u,v)  using  the 
linear  systems  of  equations.  Energy  norms  are  defined  by  ||||u||||2  m  £T(u,v).  The  discrete 
norms  |l||u||||#,  0  <  s  <  2  and  u  €  Mj,  are  defined  as  usual.  We  prove  the  following 
theorem  in  Douglas  [13]: 

Theorem  2.10:  Let  j  >  1  be  a  fixed  integer  and  d  the  dimension  of  the  problem.  Then 
the  following  holds: 

(a)  the  linear  systems  are  symmetric  and  have  Nj  positive  eigenvalues  which  are  bounded 
by  Ch]",2  where  C  is  independent  of  j. 

(b)  For  0  *-  0,  1,  or  2,  there  exist  positive  constants  Cj  ^  and  C2^  such  that 

CJIIulll^  <  llllulHI,  <  C^ltlulltfforallueM.'. 

(')  lll|Eb^en||||0  <  K(d)  C2)1  h2  ||||en||||2,  where  K(l)  -  1  and  K(2)  -  81/2. 

Theorems  1.9,  1.10,  and  1.13  can  be  rewritten  using  the  |j||-||||  norm  instead  of  the  |||  ||| 
norm.  This  proves  that  the  variable  coefficient  Dirichlet  problems  (2.4)  and  (2.5)  can  be 
solved  using  the  theory  of  Chapter  I  in  0(N]t)  operations. 
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